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Viscous-Invlscid  Interaction  in  Transonic  Flow 


Viscous- Inviscid  Interaction  in  Transonic  Flow 


by 

Laurence  B.  Wigton 

Abstract 


The  aim  of  this  thesis  is  to  couple  an  inviscid  two 
dimensional  steady  transonic  flow  calculation  with  a 
boundary  layer  calculation.  This  interaction  is  especially 
important  in  transonic  problems  since  the  boundary  layer 
has  a  significant  effect  on  the  inviscid  portion  of  the 
flow.  Here,  the  inviscid  solution  is  obtained  by  an 
algorithm  developed  for  the  full  potential  equation  by 
Holst  and  Ballhaus  while  the  attached  and  separated 
turbulent  boundary  layer  calculations  are  performed  by 
Green's  lag  entrainment  method. 

Guided  by  a  model  problem  suggested  by  Le  Balleur, 
a  viscous-inviscid  coupling  algorithm  is  developed. 
Theoretical  analysis  indicates  that  it  coverges  rapidly 
for  attached  flows  and  also  performs  well  for  separated 
flows.  These  conclusions  are  confirmed  through  a  series 
of  challenging  transonic  calculations  involving  both 
attached  and  separated  flows.  The  coupling  algorithm 
is  remarkably  stable  and  allows  computation  of  coupled 
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1.  Introduction 


The  aim  of  this  thesis  is  to  couple  an  inviscid  two 

dimensional  steady  transonic  flow  calculation  with  a 

boundary  layer  calculation.  This  interaction  is 

especially  important  in  transonic  problems  since  the 

boundary  layer  can  have  a  significant  effect  on  the 

inviscid  portion  of  the  flow.  An  example  of  a  purely 

inviscid  pressure  distribution  and  that  computed  as  the 

result  of  coupling  with  the  boundary  layer  is  shown  in 

Fig.  8.21,  As  can  be  seen,  the  change  in  the  pressure 

distribution  is  quite  dramatic.  The  viscous  effects 

smooth  out  the  pressure  near  the  trailing  edge  and 

weaken  the  shock  driving  it  upstream.  The  shock  jump 

conditions  appear  to  be  affected  as  well.  Indeed  for  the 

purely  inviscid  solution,  the  critical  pressure  coefficient 
★ 


C  shown  as  a  tick  mark  on  the  C  axis,  nearly  splits 

shock  profile.  For  the  coupled  solution,  on  the  other 

* 


the 


hand,  C  is  much  closer  to  C  downstream  of  the  shock 
P  P 

than  it  is  to  Cp  upstream  of  the  shock.  An  intuitive 
understanding  of  this  phenomenon  is  provided  by  the 
following  argument . 

For  the  purely  inviscid  case,  the  flow  must  be 
tangential  to  the  airfoil,  both  upstream  and  downstream 
of  the  shock.  Therefore,  the  shock  must  not  deflect  the 
flow,  which  forces  the  shock  to  be  normal  to  the  airfoil. 
In  the  viscous  case,  on  the  other  hand,  the  shock  subjects 
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the  boundary  layer  to  a  strongly  unfavorable  pressure 
gradient.  The  boundary  layer  responds  by  thickening  and 
forming  a  wedge  as  can  be  seen  in  Fig.  8.22.  The  wedge, 
in  turn,  forces  the  inviscid  flow  to  deflect.  In  order 
to  accommodate  the  deflecting  flow  the  shock  must  be 
oblique.  Hence  the  difference  in  the  shock  jump  conditions 
observed  in  Fig.  8.21  is  caused  by  the  change  from  a 
normal  to  an  oblique  shock. 

There  are  two  programs  in  common  use  for  computing 
coupled  viscous-inviscid  transonic  flows.  One  is  the 
viscous  Garabedian  and  Korn  (VGK)  program  developed  by 
Collyer  and  Lock  (see  Collyer  and  Lock  [1979]  and  Lock 
[1980]).  The  other  is  GRUMFOIL,  developed  by  Melnik,  Chow, 
Mead  and  Jameson  (see  Melnik  [1980]).  Both  programs  are 
capable  of  predicting  flows  around  airfoils  far  more 
accurately  than  is  possible  with  purely  inviscid  codes. 

Each  employs  the  classical  iteration  procedure  to  couple 
the  boundary  layer  with  the  outer  inviscid  flow.  This  is 
set  out  below. 


CLASSICAL  ITERATION 


Given  Geometry  Calculate 
Inviscid  Flow 


Given  Inviscid  Velocity 
Calculate  Boundary  Layer 


Use  Result  of  Boundary  Layer 
Calculation  to  Update  Geometry 
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Unfortunately,  as  is  well  known,  the  classical  iteration 
breaks  down  when  the  boundary  layer  separates  or  approaches 
separation.  (The  precise  reasons  for  this  failure  will  be 
explained  below.)  Accordingly,  the  authors  of  the  VGK 
program  and  of  GRUMFOIL  expressed  the  need  for  alternative 
coupling  procedures  to  handle  separated  flows  more  effec¬ 
tively.  Such  a  coupling  algorithm  is  developed  in  this 
thesis.  The  method  is  a  modification  of  that  proposed  by 
Le  Balleur  [1978].  The  modified  method  converges  very 
rapidly  in  regions  of  attached  flow  and  is  also  capable  of 
handling  regions  of  separated  flow.  Before  describing 
this  modification  of  Le  Balleur’ s  coupling  algorithm,  we 
will  first  discuss  the  boundary  layer  calculation  and  the 
inviscid  flow  calculation. 


2.  Green's  Lag-Entrainment  Method 


In  the  present  study  the  boundary  layer  is  calculated 
using  Green's  lag-entrainment  method  (see  Green  et  al . 

[1977]  and  East  et  al .  [1977]).  This  is  an  integral 
method  composed  of  a  coupled  system  of  first  order  ordinary 
differential  equations.  The  turbulence  modeling  is  derived 
from  the  Bradshaw-Ferriss  [1971]  turbulent  energy  equation. 
A  useful  feature  of  the  method  is  the  provision  for  per¬ 
forming  wake  calculations.  The  method  has  proven  to  be 
very  fast  and  reliable,  making  it  an  ideal  candidate  for 
use  in  a  viscous-inviscid  coupling  procedure.  The  method 

is  in  fact  used  by  the  VGK  and  GRUMFOIL  programs. 

* 

When  using  the  displacement  thickness  (6  )  concept  to 
represent  the  influence  of  the  boundary  layer  on  the  outer 
inviscid  flow,  it  is  important  to  note  that  all  integral 
boundary  layer  methods  can  be  written  in  the  form 

*  * 


A 


d6 

dx 


u  dx 


+  C 


(2.1) 


If  A  =  0,  then  it  is  impossible  to  use  (2.1)  to  solve  for 
* 

d6  /dx  in  terms  of  du/dx.  In  this  case  it  is  not  appropri¬ 
ate  to  integrate  the  boundary  layer  equations  using  a 
prescribed  value  of  du/dx.  This  is  referred  to  as  the 
Goldstein  or  separation  singularity.  The  separation 

singularity  can  be  avoided  by  integrating  the  boundary 

* 

layer  equations  with  a  prescribed  value  of  d6  /dx  and 
solving  for  du/dx.  This  is  the  so-called  inverse  mode  of 


5 


calculating  the  boundary  layer  flow.  On  the  other  hand, 

if  B  =  0  somewhere  in  the  interval  of  integration,  then  it 

becomes  impossible  to  use  (2.1)  to  solve  for  du/dx  in 
* 

terms  of  d5  /dx.  In  this  case  it  becomes  inappropriate 

to  integrate  the  boundary  layer  equations  with  prescribed 
* 

values  of  d6  /dx.  This  situation  is  referred  to  as  the 
Crocco-Lees  singularity. 

A  sharper  criterion  for  deciding  whether  to  prescribe 

♦ 

d6  /dx  or  du/dx  along  the  boundary  layer  edge  is  provided 
by  considerations  of  numerical  accuracy.  If  one  is  inte¬ 
grating  numerically  a  differential  equation  of  the  form 

^  =  Dy  (2.2) 


then  the  accuracy  of  the  result  depends  on  the  size  of  Dh 
where  h  is  the  step  size.  To  ensure  numerical  accuracy 
we  would  like  to  'arrange  for  D  to  be  as  small  as  possible. 
Equation  (2.1)  can  be  rewritten  in  either  of  the  forms 

* 


d6  =  /B  i  4.  £ 

dx  'A  u  dx'^  A 


(2.3) 


du  =  /A  _1^  d6*..  C 

dx  '•B  dx  ■  B 

0 


(2.4) 


* 

Comparing  the  coefficient  of  6  in  Eq.  (2.3)  with  the 
coefficient  of  u  in  Eq.  (2.4)  we  see  that  from  the  view¬ 
point  of  numerical  accuracy,  it  is  better  to  prescribe 
du/dx  at  the  boundary  layer  edge  when 

lA  d6*  ■  iB  1^  dui 

' B  ^  *  dx  ^  ' A  u  dx  ^ 


(2.5) 
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* 

otherwise  it  is  better  to  prescribe  d6  /dx. 

It  is  shown  by  East  [1977]  that  Green's  lag-entrainment 
equations  imply  that 


dS  =  §.^=17  I  ^2  6  du 

dx  1  ^2  u  dx  ^1  H  u  dx 


(2.6) 


where  6  is  the  momentum  thickness.  The  third  expression 

* 

in  Eq.  (2.6)  merely  serves  to  define  H=6  /6.  A  plot  of 
F2/H  versus  a  shape  parameter  H  for  both  the  standard 
Green's  lag-entrainment  method  and  for  a  modified  method 
due  to  East  is  shown  in  Fig.  2.1.  Separation  normally 
occurs  for  H  >  2.4  (approximately).  For  H  <  1.6  the 
modified  and  standard  method  are  identical.  The  modified 
method  was  introduced  by  East  to  improve  the  results 
based  on  Green's  method  for  separated  flows.  As  can  be 
seen  from  Fig.  2.1,  for  separated  flows  this  modification 
makes  the  value  of  IF2/HI  larger  than  its  corresponding 
value  yielded  by  the  standard  method.  For  Mach  numbers 
in  the  range  of  transonic  interest,  Fg/H  is  never  0  or  » . 
In  fact  1’2/H  is  always  a  well-defined  negative  number. 

The  Crocco-Lees  singularity  which  corresponds  to  ^2/^  =  0 
does  not  occur  until  the  Mach  number  exceeds  1.5.  The 


separation  singularity  which  corresponds  to  F2/H  =  >»  never 
occurs.  However,  in  view  of  the  analysis  given  above,  if 
we  are  to  expect  accurate  results  when  integrating  the 
boundary  layer  equations  with  a  prescribed  velocity 
gradient ,  we  must  have 


I2  1  ^ 

H  u  dx 


hi 


<<  1 


(2.7) 
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where  h  is  again  the  step  size.  This  criterion  is  dif¬ 
ficult  to  meet  for  separated  flows  because  IF2/HI  attains 
such  large  values.  This  is  all  the  more  true  for  shock 
induced  separation  because  then  Idu/dxj  is  also  large. 
Since  the  classical  coupling  procedure  depends  on  inte¬ 
grating  the  boundary  layer  equations  with  a  prescribed 
value  of  du/dx,  the  inequality  (2.7)  explains  why  trouble 
can  be  expected  for  separated  flows  even  though  there  is 
no  separation  singularity.  This  demonstrates  the  need  to 
develop  coupling  procedures  which  allow  one  to  integrate 
the  boundary  layer  equations  with  a  prescribed  d6  /dx. 

An  additional  reason  for  the  failure  of  the  classical 
iteration  will  be  given  when  we  discuss  Le  Balleur's 
method. 


Using  the  same  reasoning  which  led  inequality 

(2.5)  we  see  that  we  would  prefer  to  integrate  Green's 

♦ 

boundary  layer  equations  with  a  prescribed  d6  /dx  whenever 


iJL  A.  ^§1\  ^  i  ^1 

'F2  dx  ‘  '  H  u  dx' 


(2.8) 


For  typical  airfoil  calculations  Eq.  (2.8)  usually 

indicates  that  Green’s  lag-entrainment  equations  should  be 

* 

integrated  with  a  prescribed  value  for  d6  /dx  everywhere 
except  in  the  far  wake  region. 

In  both  the  VGK  and  GRUMFOIL  programs  the  effect  of 
the  boundary  layer  on  the  inviscid  flow  is  represented 
through  a  surface  transpiration  condition 


(2.9) 
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In  this  formula 

V  =  normal  velocity 
n 

=  density 

=  tangential  velocity 

^  =  derivative  in  tangential  direction 
ds 

According  to  the  analysis  presented  by  Le  Balleur 

[1977],  use  of  Eq.  (2.9)  avoids  the  Crocco-Lees  singularity; 

* 

that  is,  if  we  define  m = p^u^S  and  use  the  integral 
boundary  layer  equations  to  derive  an  equation  of  the 
form 

+  (2.10) 

dx  u  dx 

then  the  coefficient  B'  never  vanishes  (regardless  of 
Mach  number).  This  fact  is  easily  verified  for  Green’s 
lag-entrainment  method  by  numerical  computation. 

Another  advantage  of  the  surface  transpiration 
condition  (2.9)  is  that  it  does  not  require  the  regenera¬ 
tion  of  a  mesh  after  each  iteration,  as  would  be  necessary 
if  the  displacement  thickness  concept  were  implemented 
exactly . 

In  the  present  study,  both  the  displacement  thickness 
and  the  surface  transpiration  concepts  are  implemented. 
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3.0  Inviscid  Analysis 

Calculation  of  steady  inviscid  transonic  flows  is  a 
formidable  task.  The  equations  governing  transonic  flow 
are  inherently  nonlinear  and  change  type  within  the  solu¬ 
tion  domain,  from  elliptic  in  subsonic  regions  to  hyper¬ 
bolic  in  supersonic  regions.  Moreover,  one  must  provide 
for  embedded  shocks,  the  positions  of  which  are  to  be 
calculated  as  part  of  the  solution.  The  most  successful 
numerical  methods  for  calculating  inviscid  transonic  flow 
on  a  routine  basis  use  finite  difference  schemes.  There¬ 
fore,  only  these  methods  were  considered  for  use  in  the 
viscous-inviscid  coupling  problem. 

Magnus  and  Yoshihara  [1970]  solved  the  unsteady  Euler 
equations  using  an  explicit  second-order  difference  scheme 
similar  to  that  developed  by  Lax  and  Wendroff  (see  Lax 
[1954]).  The  unsteady  equations  are  always  hyperbolic. 
Thus  the  problem  of  dealing  with  equations  of  mixed  type 
was  circumvented.  The  difference  scheme  effectively  added 
a  numerical  viscosity  which  allowed  shocks  to  appear  as 
rapid  but  continuous  changes  in  the  flow.  This  allowed 
shocks  to  be  captured  as  a  natural  part  of  the  solution. 
The  disadvantages  of  the  method  of  Magnus  and  Yoshihara 
are  two-fold.  Firstly,  the  equations  were  solved  in 
Cartesian  coordinates  which  forced  them  to  employ  a 
cumbersome  set  of  embedded  meshes  to  enforce  the  boundary 
conditions  on  the  airfoil.  Secondly,  the  rate  at  which 
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the  unsteady  solution  approaches  the  desired  steady  state 
solution  is  very  slow.  Therefore,  the  method  is  computa¬ 
tionally  expensive.  Of  course,  considerable  progress  has 
since  been  made  in  solving  the  Euler  equations  (see  for 
example  Beam  and  Warming  [1976]).  However  the  calcula¬ 
tions  are  still  very  expensive.  Therefore  it  was 
decided  not  to  solve  the  Euler  equations  as  part  of  our 
viscous-inviscid  coupling  problem,  but  to  rely  on  the  full 
potential  equations. 

A  significant  forerunner  to  modern  methods  for  solving 
the  latter  was  provided  by  Murman  and  Cole  [1971].  In 
this  the  small  disturbance  eq_ations  were  solved  with  cor¬ 
responding  boundary  conditions.  Type  dependent  differencing 
was  introduced.  Centered  differences  were  used  in  elliptic 
regions  while  upwind  differences  were  employed  in  hyper¬ 
bolic  regions.  A  special  shock  point  operator  was  incor¬ 
porated  to  preserve  the  conservative  nature  of  the  scheme. 
The  difference  equations  were  solved  by  successive  line 
over-relaxation  (SLXDR).  The  resulting  calculation  is 
faster  than  the  methods  for  solving  the  Euler  equations 
by  several  orders  of  magnitude. 

Encouraged  by  the  success  of  Murman  and  Cole,  Jameson 
[1974]  developed  methods  for  solving  the  non-conservative 
form  of  the  full  potential  equation: 

-  0  (3.1) 

(see  for  example  Shapiro  [1953]).  A  circle  plane  mapping 
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which  transformed  the  exterior  of  the  airfoil  to  the 
interior  of  the  unit  circle  was  employed.  This  allowed 
the  exact  boundary  conditions  at  the  airfoil  to  be  easily 
enforced.  Upwind  differencing  in  the  direction  of  the 
flow  was  used.  Also  artificial  time  terms  guaranteed  the 
stability  of  the  method.  The  difference  equations  were 
also  solved  using  SLOR. 

} 

Later  Jameson  [1975]  extended  his  method  to  the 
conservative  form  of  the  full  potential  equations: 

(p$  )  + (p$  )  =  0  (3.2) 

This  conservative  form  is  to  be  preferred  over  the  non¬ 
conservative  form  because  it  assures  the  proper  jump  con¬ 
ditions  across  shocks.  If  the  non-conservative  form  (3.1) 
is  used,  then  the  shock  location  is  mesh  dependent. 

In  an  attempt  to  improve  the  rate  of  convergence  to 
the  solution  of  inviscid  transonic  flows,  Ballhaus, 

Jameson  and  Albert  [1978]  developed  an  implicit  approxi¬ 
mate  factorization  (AF)  algorithm.  The  method  was  applied 
to  the  steady  state  transonic  small  disturbance  equation 
and  enjoyed  a  rate  of  convergence  5-7  times  as  great  as 
that  achievable  by  the  SLOR  algorithm.  Holst  and  Ballhaus 
[1979]  followed  up  on  this  success  by  applying  the  (AF) 
algorithm  to  the  conservative  full  potential  equation  in 
Cartesian  coordinates.  Holst  [1979]  later  extended  the 
f  method  to  a  numerically  generated  body  fitted  coordinate 

system,  which  allows  the  exact  boundary  conditions  at  the 


f 
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airfoil  to  be  implemented  easily. 

The  fastest  method  for  solving  the  full  potential 
equation  numerically  was  developed  by  Jameson  [1979] .  In 
this  a  multiple  grid  method  developed  by  Achi  Brandt  [1979] 
is  combined  with  an  Alternating  Directions  Implicit  (ADI) 
scheme  similar  to  the  (AF)  scheme  developed  by  Holst  and 
Ballhaus.  The  ADI  scheme  acts  as  the  smoothing  algorithm 
necessary  for  the  success  of  the  Multiple  Grid  method.  The 
resulting  Multigrid  Alternating  Direction  (MAD)  method  con¬ 
verges  very  rapidly  indeed.  A  typical  calculation  requires 
only  about  10  iterations  to  achieve  convergence. 

For  the  purpose  of  calculating  a  coupled  viscous- 
in  vise  id  flow  it  was  decided  not  to  use  Jameson's  MAD 
method.  This  decision  was  based  on  the  difficulty  in 
implementing  multigrid  methods.  More  importantly,  it  is 
difficult  to  construct  a  coupling  algorithm  which  is  suf¬ 
ficiently  rapid  to  take  advantage  of  the  speed  with  which 
MAD  calculates  the  inviscid  flow.  Thus  using  MAD  would  not 
necessarily  improve  the  rate  at  which  a  coupled  viscous- 
inviscid  calculation  can  be  carried  out.  Accordingly,  it 
was  decided  to  use  the  (AF)  scheme  developed  by  Holst. 

This  method  is  also  very  fast  and  reliable.  Moreover,  it 
is  well  documented  allowing  it  to  be  easily  implemented  on 
the  computer. 

3. 1  Motivation  for  Approximate  Factorization  Schemes 

An  intuitive  motivation  for  approximate  factorization 
(AF)  schemes  is  provided  by  the  following  considerations. 


13 


Suppose  that  it  is  required  to  solve  a  problem  of 
the  form 


L<j!  =  0  (3.3) 

where  L  is  an  operator.  Given  an  approximation  <()”  to  <p 
and  an  approximation  N  to  L,  we  might  be  able  to  find  a 
better  approximation  to  <p  by  solving 

-()>”)  =  -  (3.4) 

Indeed,  in  the  case  where  L  is  linear  and  N  =  L,  then 

y.  I  T 

solving  Eq.  (3.4)  for  yields  the  exact  solution  to 

Eq.  (3.3).  We  can  introduce  over-relaxation  into  the 
algorithm  indicated  by  Eq.  (3.4)  by  solving  for 

..n+l.,*  _  /-o  CN 

(<J)  )  -  oj<p  +  (l-ai)ip  (3.5) 

which  implies  that 

(cf)  )-<()  =(0((j>  -<)>)  (3.6) 

If  N  is  linear,  then  Eqs.  (3.4)  and  (3.6)  yield 

N(((f"'^^)’*‘ -  -  ())")  =-a)L4>"  (3.7) 

Thus  the  version  of  Eq.  (3.4)  which  includes  over- relaxa¬ 
tion  can  be  written 

N(<!)'^'^^  -  4)")  =  -coLfl)"  (3.8) 

In  this  formula,  is  the  residual  which  is  a  measure  of 

"t 

how  well  the  operator  is  satisfied  by  the  n^  level  solution 
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I 

(l”,  oj  is  the  relaxation  parameter  and  is  the 

correction.  The  iteration  given  by  Eq .  (3.8)  can  be 
regarded  as  an  iteration  in  pseudo-time,  where  the  super¬ 
script  n  indicates  the  time  level  of  the  iteration.  In 
the  approximate  factorization  scheme,  N  is  chosen  to  be 
the  product  of  two  or  more  factors  indicated  by 

N  =  (3.9) 

The  factors  and  N2  are  chosen  so  that  they  are  easily 
invertible,  their  product  is  an  approximation  to  L,  and 
the  iteration  indicated  by  Eq .  (3.8)  is  stable. 

3.2  Artificial  Density 

When  solving  the  full  potential  equations  in  conserva¬ 
tion  form  Eq.  (3.2),  the  shock  jump  conditions  are 

fp$  In  +  [p$  In  =0  (3.10) 

x^  X  y-*  y  ^  ^ 

) 

I 

!  In  this  formula  (n^.n^)  is  normal  to  the  shock  and  [  ] 

I  denotes  the  change  in  a  quantity  when  crossing  the  shock. 

Unfortunately,  Eq .  (3.10)  is  not  enough  to  ensure  unique¬ 
ness  of  solutions.  Non-admissible  expansion  shocks,  for 
example,  are  not  excluded  by  Eq.  (3.10). 

One  'nethod  to  ensure  uniqueness  is  to  add  viscous 
terms.  iscous  terms  occur  naturally  in  the  Navier-Stokes 
equations.  They  guarantee  that  shocks  are  regions  of  rapid 
but  continuous  change  and  exclude  the  possibility  of  expan¬ 
sion  shocks.  When  solving  the  inviscid  flow  equations 
involving  shocks,  it  is  customary  to  add  numerically 
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convenient  viscous  terns  which  smear  shocks  over  three  to 
four  mesh  points  and  also  exclude  physically  unrealistic 
solutions. 

In  earlier  attempts  to  solve  the  inviscid  transonic 
flow  equations,  artificial  viscosity  was  explicitly  added 
through  type  dependent  difference  schemes.  Centered  dif¬ 
ferences  were  used  in  elliptic  regions  and  upwind  dif¬ 
ferences  were  used  in  hyperbolic  regions.  Recently  it  has 
been  discovered  that  it  is  more  convenient  to  add  artificial 
viscosity  by  retarding  the  density.  The  method  used  in  the 
present  study  is  to  use  p  given  by 

p  =  p  _  y  (3.11) 

3  s 


where : 


=  derivative  of  p  in  the  streamwise  direction  (3.12) 

0  S 


As  =  mesh  width 


(3.13) 


y  = 


= 


Min[C(.\r-l),l]  M  >  1 

0  M  <  1 


(3.14) 


M  =  Mach  number 


C  =  bias  coefficient,  usually  1.5 


(3.15) 

(3.16) 


For  completeness  we  mention  that 


v2  *2  ^  .2 

V  =  + 


=  (l-v2)l/<^-l) 


(3.17) 


(3.18) 
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V^/(l-V^)  (3.19) 

Y  =  ratio  of  specific  heats  =  1.4  for  air  (3.20) 

where  the  velocity  V  is  normalized  by  the  maximum  velocity 
and  the  density  p  is  normalized  by  the  stagnation 
density . 

Numerical  experiments  reveal  that  in  order  to  main¬ 
tain  stability,  3p/3s  must  be  calculated  using  upwind 
differencing.  Also  pressure  overshoots  upstream  of  the 
shock  can  be  avoided  by  using  the  Mach  number  M  computed 
one  mesh  point  upstream  in  Eq.  (3.14). 

From  the  above  formulas  it  is  evident  that 

p  =  p  when  M  <  1  (3.21) 

while  for  supersonic  flow  p  is  retarded  in  the  upstream 
direction  by  an  amount  which  increases  with  Mach  number. 
Formula  (3.14)  ensures  that  p  is  never  retarded  for  more 
than  one  mesh  point.  In  place  of  Eq.  (3.2)  we  now  solve 

(P'l'  )  +  (p^  )  =  0  (3.22) 

^  x'x  ^  y  y 

The  new  jump  conditions  become 

tP^ylOy  =  0  (3.23) 

From  Eq .  (3.11)  we  see  that  as  As->-0,  p -*■  p .  Thus  by 
comparing  Eq.  (3.23)  with  Eq.  (3.10)  we  see  that  the  exact 
shock  jump  conditions  are  enforced  as  the  mesh  width  tends 
to  0. 
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The  artificial  density  method  is  a  convenient  way  of 
introducing  artificial  viscosity  which  preserves  the 
proper  shock  jump  conditions.  Moreover,  as  we  shall  see, 
it  permits  the  full  potential  equations  to  be  solved 
using  an  approximate  factorization  scheme. 


3 . 3  Approximate  Factorization  Applied  to  the  Full 
Potential  Equation 

In  order  to  solve  Eq.  (3.22)  we  introduce  the  finite 
difference  analog: 


L({i  =  (6  P  •  .  1  ■  &  +  <S  p.  -.1  5  )4> . 


(3.24) 


The  first  factorization  introduced  by  Ballhaus,  Jameson 
and  Albert  [1978]  referred  to  as  AFl  corresponds  to 

•4-  -►  *4-  -♦» 

aN  =  -  (a  -  5^  p^^^  -  6^  p  6^)  (3.25) 


where  a  is  a  positive  number  v/hich  is  part  of  a  parameter 
sequence  chosen  so  as  to  optimize  the  rate  of  convergence 
of  the  algorithm  (3.8).  In  general  large  values  of  a 
reduce  high  frequency  errors  while  small  values  of  a 
reduce  low  frequency  errors.  The  best  sequence  of 
parameters  to  use  for  any  particular  problem  must  be 
determined  experimentally  since  the  problem  is  non-linear. 
However,  Ballhaus  et  al .  [1978]  do  provide  rough  guide¬ 

lines  based  on  a  linear  analysis.  Fortunately,  almost  any 
reasonable  choice  of  parameter  si  ^  ence  leads  to  an  algo¬ 
rithm  which  converges  much  faster  than  SLOR. 

The  AFl  scheme  Indicated  in  Eq .  (3.25)  performs  very 
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well  for  subsonic  flows  but  not  for  supersonic  flows.  A 
heuristic  argument  explaining  this  behavior  will  now  be 
presented . 

If  we  regard  to  be  (j>^,  then  the  left-hand 

side  of  Eq .  (3.8)  with  N  given  by  Eq .  (3.25)  contains  a 
term.  In  regions  of  subsonic  flow  the  right-hand  side 
of  Eq.  (3.8)  is  an  elliptic  operator,  whereas  in  regions 
of  supersonic  flow  it  is  hyperbolic.  In  subsonic  regions 
we  can  model  the  AFl  algorithm  with  the  equation 


<f> 


t 


XX 


yy 


(3.26) 


while  for  supersonic  flow  we  use  the  model  equation 


t 


Assume  solutions  of  the  form 


(3.27) 


=  g(7t+i(ax+by)  (3 

where  a  and  b  are  real  constants.  Substitute  Eq.  (3.28) 
into  Eqs.  (3.26)  and  (3.27)  and  solve  for  a.  We  find 


0  = 


2 


) 


for  subsonic  flow 
for  supersonic  flow 


(3.29) 


In  'the  subsonic  case  we  see  that  (})  always  remains  bounded 
with  time  but  this  is  not  true  in  the  supersonic  case. 
This  explains  the  behavior  of  the  AFl  scheme.  The  <})^ 
term,  which  is  implicitly  introduced,  is  compatible  with 
subsonic  flows  but  not  with  supersonic  flows. 
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Time  dependent  terms  of  the  form  4)^^  or  on  the 
other  hand,  are  compatible  with  both  subsonic  and  super¬ 
sonic  flows.  Indeed,  if  we  consider 


xt 

=  4)  +  4> 

XX  yy 

(3.30) 

xt 

=  4>  -  4> 

^xx  ^yy 

(3.31) 

and  assume  solutions  of  the  form  given  in  Eq.  (3.28),  then, 
in  both  cases,  we  find  that  a  is  purely  imaginary  so  that  <j> 
is  bounded  with  time. 

Accordingly,  Ballhaus  et  al .  [1978]  developed  a  second 
approximate  factorization  scheme  called  AF2  which  intro¬ 
duces  the  proper  time  dependent  terms.  The  version  of  AF2 
used  in  the  present  study  is  given  by 

aN  «  (a4y  Pj+i)(a5y4^  (3.32) 

In  order  to  enhance  further  the  stability  of  AF2 ,  a  term 
of  the  form 

±  a  8  I  (3.33) 

is  added  to  the  second  factor  for  N  given  in  Eq .  (3.32). 

For  subsonic  regions  6  is  usually  assigned  the  value  0.3. 

For  supersonic  regions  6  is  a  user  defined  constant 
attaining  values  as  large  as  9.0  for  particularly  chal¬ 
lenging  flows.  The  double  arrow  notation  indicates  that 
the  difference  is  always  upwind  and  the  sign  is  chosen  so 
as  to  increase  the  magnitude  of  the  diagonal  terms  in  the 
matrix  for  the  second  factor  of  N.  The  first  factor  for  N 
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(3.37) 
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—  ^  e  / —  k  9,  9  ^  j 

-  J  33^  (P  /i  g  53^)w  dy 

_  13^  f —  kJl  9$  ^  , 

Comparing  the  first  and  last  expressions  we  conclude  that 
in  y  coordinates,  the  full  potential  equation  is  given  by 

9  (p  /g  gk.  ^  0  (3  4,^ 

iL  K 

Once  a  mesh  has  been  generated,  the  x  coordinates  are 


known  in  terms  of  the  y  coordinates.  This  allows  us  to 
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compute  from  Eq .  (3.37)  and  hence  (g^*^)  from  Eq. 

(3.38)  and  g  from  Eq.  (3.39).  This  provides  us  with  the 
metric  coefficients  required  in  Eq.  (3.40).  Also  in  y 
coordinates  we  have 


„2  _  .  3<I>  3$  _  „ij  3$  3$ 

-  «ij  3ir  337  -  e  ^  ^ 


(3.41) 


allowing  us  to  compute  p  and  M  . 

Naturally,  in  solving  Eq.  (3.40)  we  replace  p  with  p 
and  use  the  finite  difference  analog  of  Eq.  (3.40).  The 
precise  algebraic  details  may  be  found  in  Holst  [1979] . 
The  AF2  scheme  used  in  the  present  study  was  the  finite 
difference  version  of 


3l;*3|7  <3-«) 

The  finite  difference  operators  were  chosen  so  as  to  make 
the  diagonal  terms  as  large  as  possible  in  absolute  value. 

Two  versions  of  the  code  exist.  The  first  uses 
Cartesian  coordinates  and  applies  the  boundary  conditions 
according  to 

^y/$x  =  ^  y  =  °  (3.43) 

where  f  is  the  slope  of  the  airfoil  and  d6  /dx  is  the 
slope  of  the  boundary  layer  displacement  thickness.  This 
version  of  the  code  is  intended  for  testing  various  viscous- 
inviscid  coupling  techniques. 

The  more  exact  version  of  the  code  uses  Eisman's 
[1979]  method  to  generate  a  body  fitted  mesh  and  imposes 
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the  boundary  condition 


13  * 

V  =  —  ^  (p  u  S  ) 
n  p  3s  e  e 
e 


(3.44) 


at  the  airfoil  surface  and  along  the  wake  centerline. 
Eisman's  method  was  chosen  to  generate  the  mesh  because 
it  is  very  fast  and  gives  excellent  control.  A  typical 
mesh  is  shown  in  Figs.  3.1  and  3.2.  Note  that  the  mesh  is 
designed  to  conform  with  both  the  symmetric  airfoil  and 
the  wake  at  the  centerline.  This  makes  it  easy  to 
implement  the  boundary  conditions  ne.cessary  to  solve  the 
viscous-inviscid  coupling  problem. 
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4.  Le  Balleur's  Model  Problem 

Although  Le  Balleur's  [1978]  method  for  analysing 
viscous-inviscid  coupling  algorithms  is  relatively  easy,  it 
has  not  yet  received  the  recognition  in  this  country  that 
it  deserves.  Therefore,  some  aspects  of  his  theory  will  be 
presented  in  detail. 

In  order  to  gain  insight  into  coupling  algorithms, 

Le  Balleur  considers  a  model  problem.  The  inviscid  flow 
is  assumed  to  be  governed,  over  the  half  plane  y>0,  by 
the  linearized  small  disturbance  equation 

(l-M^)<^^j^+4>yy  =  0  ,  (4.1) 

where  the  Mach  number  M  is  constant  and  the  full  potential 
$  is  given  by 

$  =  u^(x  +  (i))  (4.2) 

so  that  the  inviscid  velocity  u  is  given  by 


u  =  =  u^(l  +  0^) 


dx  «  XX 


(4.3) 

(4.4) 


The  boundary  layer  is  assumed  to  be  governed  by  the  equation 


^  =  A  +  B  ^ 
dx  dx 


(4.5) 


where  A  and  B  are  constants.  For  convenience  we  introduce 

* 


a  = 


d6 

dx 


(4.6) 


i 


I 
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and  rewrite  the  boundary  layer  equation  as 

a  =  A  +  B  (4.7) 

The  boundary  conditions  for  the  inviscid  flow  are  applied 
in  a  small  disturbance  fashion  so  that 

(t)y=f'+Qi  at  y=0  (4.8) 

where  f  is  the  slope  of  the  body  about  which  we  are  com¬ 
puting  the  flow. 

We  will  have  a  coupled  viscous-inviscid  solution  if  we 
can  determine  a  so  that  the  inviscid  velocity  gradient 
(evaluated  at  y  =  0)  agrees  with  the  velocity  gradient  com¬ 
puted  from  the  boundary  layer  equation.  Suppose  that  such 
a  distribution  of  a  is  given  to  us.  In  order  to  analyse  a 
particular  algorithm  for  calculating  a  coupled  solution, 

i.  V  X 

Le  Balleur  perturbs  the  exact  a  by  Aa  =  ee  .  He  then 
applies  the  coupling  algorithm  to  see  if  the  error  in  a  is 
increased  or  decreased.  The  analysis  of  the  classical 
iteration,  for  example,  is  now  described. 

If  we  perturb  the  exact  a  by  an  amount  Aa  =  ee^^’^,  it 
follows  from  Eqs.  (4.1)  and  (4.8)  that  the  exact  ()>  will  be 
perturbed  by  A({)  governed  by 

(l-M2)(A<t)^j^  +  (A4))yy  =  0  (4.9) 

where 

(A()))y  =  ee^^'^  at  y  =  0  (4.10) 

Introduce 


I 
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B  =  1-M 


(4.11) 


If  the  inviscid  flow  is  subsonic,  then  Eq.  (4.9)  may  be 


rewritten  as 


B  (A(j))  +  (All))  =  0 

XX  ^  ^yy 


(4.12) 


The  solution  to  Eq.  (4.12),  subject  to  the  boundary  con¬ 
dition  (4.10),  is 


All  = 


^^ivx^-Bvy 


(4.13) 


The  solution  involving  is  discarded  because  it  is  not 

well  behaved  for  y -»■  «> .  Now,  using  Eq.  (4.4),  we  see  that 
this  perturbation  of  41  causes  the  inviscid  velocity  gradient 
(at  y = 0)  to  be  changed  by  the  amount 


j  u  V  . 

^(dJ>  -  — 


(4.14) 


The  velocity  gradient  which  the  inviscid  flow  imposes  on  the 
boundary  layer  differs  from  the  exact  velocity  gradient  by 
the  amount  shown  in  Eq.  (4.14).  Thus,  by  Eq.  (4.7),  the 
boundary  layer  is  forced  to  respond  by  giving  a  value  of  a 
which  is  in  error  by 


vBu 

,  ”  ivx 

Aa  =  — g —  ee 


(4.15) 


Thus  in  the  subsonic  case,  the  classical  iteration  mutliplies 
the  original  error  in  a,  namely  Aa  =  ee^'^^,  by  a  factor 


for  M  <  1 


(4.16) 


Similarly  for  the  supersonic  case  we  find 
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-0^(A(j))  +  (A<p)  =0  (4.17) 

XX  yy 

ivx  -iBvy 

AO,  =  ie - e - ^ 

-iBv 

Only  the  right  running  wave  is  considered  in  Eq .  (4.18). 
Following  the  same  steps  as  those  used  to  derive  Eq .  (4.16) 
we  find 

-ivBu 

M  =  - 5—^  for  M  >  1  (4.19) 

p 

The  classical  iteration  will  converge  only  if  1m1  <  1 
for  all  possible  values  of  v.  The  latter  lie  in  the 
approximate  range 


i  <  V  < 


Ax 


(4.20) 


where  Ax  is  the  mesh  width  and  L  is  the  size  of  the  computa¬ 
tional  domain.  Also  we  should  note  that  for  separated  flows, 
Bu^  becomes  large.  Thus,  for  separated  flows  computed  on 
fine  grids,  we  can  easily  expect  that  |ul  >  1  in  either  the 
subsonic  or  supersonic  case.  This  provides  the  second 
reason  for  expecting  the  classical  iteration  to  fail  for 
separated  flows.  If  we  refine  the  mesh  sufficiently  so  that 
we  can  accurately  integrate  the  boundary  layer  equations 
with  a  prescribed  velocity  gradient  when  the  flow  separates, 
we  can  easily  force  | y |  >  1 ,  so  that  the  viscous-inviscid 
coupling  will  diverge. 

It  is  of  interest  to  see  how  the  classical  iteration 
operates  when  under- relaxation  is  employed.  For  this 
purpose  it  is  useful  to  note  that,  according  to  the  above 
analysis,  the  classical  iteration  can  be  viewed  as  a 
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linear  operator  £  with  eigenfunctions  e 


ivx 


IVXv  ,  IVX 

jC  (e  )  =  ue 


(4.21) 


where  the  eigenvalue  p  is  given  by 

vBu 


h  =  < 


-ivBu_ 


6 


for  M  <  1 


for  M  >  1 


(4.22) 


When  relaxation  is  employed  with  relaxation  factor  oj  the 
linear  operator  £  is  replaced  by 


=  (oJC  +  (l-w)I 


(4.23) 


where  I  is  the  identity  operator.  Note 


^0  =  ' 


(4.24) 

(4.25) 


From  (4.21)  and  (4.23)  we  see  that 


£  (e^^^)  =  [up  +  ( 1-00 )  ]  e 


ivx 


(4.26) 


ivx 

Thus  e  is  also  an  eigenfunction  of  £  with  corresponding 

u 

eigenvalue 


p^  =  up  +  ( 1-u  ) 


(4.27) 


Note  that 


Pq  =  1 


Pj^  =  p 


(4.28) 

(4.29) 
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As  (i)  decreases  from  1  to  0,  y  varies  from  y  to  1 . 

(i) 

Geometrically,  if  0  <  oj  <  1 ,  then  the  relationship  between 
y,  y  and  1  is  as  shown  in  Fig.  4.1.  We  can  guarantee 

CO 

the  convergence  of  the  classical  iteration  i f  we  can  choose 
OJ  so  that  a(£  ),  the  set  of  eigenvalues  of  £  ,  lies  within 

0)  CO 

the  unit  circle.  From  (4.22)  we  see  that  if  M>  1  then  con¬ 
ditions  are  as  shown  in  Fig.  4.2.  It  is  geometrically 
clear  that  if  '-e  make  co  sufficiently  small,  we  can  arrange 
for  a(£  )  to  lie  within  the  unit  circle  and  thus  ensure 
the  convergence  of  the  classical  iteration  in  this  case. 
Similar  remarks  apply  when  M>1  and  B>0.  When  the  flow  is 
subsonic,  conditions  are  as  shown  in  one  or  other  of  Figs. 

4 . 3  and  4.4. 

In  the  first  case  we  can  choose  co  so  that  a(f  )  lies 

CO 

within  the  unit  circle,  but  in  the  second  case  th’.s  is  not 
possible . 

Summarizing  these  results  we  see  that  for  supersonic 
flow  the  classical  iteration  can  always  be  stabilized  by 
under-relaxation.  For  subsonic  flow  the  classical  itera¬ 
tion  can  be  stabilized  by  under-relaxation  when  B  <  0  but 
not  necessarily  if  B  >  0. 

For  Green's  lag-entrainment  method  we  always  have 
B  <  0,  so  that  the  classical  iteration  can  always  be  stabi¬ 
lized  by  under-relaxation.  However,  it  is  still  not  a  good 
idea  to  use  the  classical  iteration  when  the  flow  separates 
because  in  this  case  we  cannot  integrate  the  boundary  layer 
equations  accurately  with  a  prescribed  velocity  gradient. 
Moreover,  the  under-relaxation  parameter  can  become  so 
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small  that  this  method  is  rendered  too  slow  for  practical 
use.  We  really  must  search  for  a  coupling  algorithm  which 

allows  us  to  integrate  the  boundary  layer  equations  with  a 

* 

prescribed  value  of  d6  /dx. 

Perhaps  the  most  natural  coupling  algorithm  which 

allows  one  to  integrate  the  boundary  layer  equations  with 

♦ 

a  prescribed  value  of  do  /dx  is  the  so-called  inverse 
method : 

Inverse  Iterat ion 


Since  we  can  associate  a  linear  operator  jC  with  the 
classical  iteration,  we  can  associate  the  inverse  operator 
with  the  inverse  iteration.  If  we  apply  X  ^  to  both 
sides  of  Eq.(4.21).,  we  find  that 


•1  /  i  vx , 

(e  )  = 


Mie 


lUX 


(4.30) 


where 
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/ 

3 

OO 


13 


for 


for 


M  <  1 


M  >  1 


(4.31) 


Thus  has  the  same  eigenfunctions  as  £  with  corresponding 

eigenvalues  which  are  the  reciprocals  of  those  of  £.  The 
largest  values  for  IPjI  occur  for  low  frequencies.  In  view 
of  Eqs.  (4.20)  and  (4.31)  it  follows  that 


=  TbT^  ^ 

Thus  max  lvij|  is  proportional  to  the  size  of  the  computa¬ 
tional  domain  L.  In  cases  where  max  luj]  >>  1,  the  size  of 
the  relaxation  parameter  w  required  to  produce  convergence 
should  vary  inversely  with  L. 

Confirmation  of  this  result  is  provided  by  a  calcula¬ 
tion  performed  by  Melnik  [1976] .  In  this  the  boundary 

layer  equations  are  integrated  using  a  prescribed  value 
* 

of  d(S  /dx  to  yield  du/dx.  Since  the  flow  considered  by 
Melnik  is  incompressible,  the  inviscid  flow  can  be  solved 
with  a  prescribed  distribution  of  du/dx  using  a  Hilbert 
integral.  Convergence  of  this  inverse  iteration  scheme  was 
first  achieved  over  a  computational  domain  [0,3]  using 
0)  =  .15.  When  the  size  of  the  computational  domain  was 
increased  to  [0,20.791]  Melnik  found  it  necessary  to 
reduce  oj  to  the  value  .02.  Hence  a  seven-fold  increase  in 
the  size  of  the  computational  domain  did  indeed  require 
roughly  a  seven-fold  decrease  in  the  size  of  co . 


32 


In  general  the  size  of  the  relaxation  parameter 
required  to  stabilize  the  inverse  iteration  over  a  reason¬ 
able  sized  computation  domain  is  disappointingly  small. 
Also  the  inverse  iteration  requires  that  the  inviscid  flow 
be  computed  with  a  prescribed  value  of  du/dx.  This 
requires  performing  a  design  calculation  which  is 
particularly  difficult  in  transonic  flow  involving  shock 
waves.  Accordingly,  Le  Balleur  develops  a  semi-inverse 

method  which  allows  one  to  calculate  the  boundary  layer 

♦ 

with  a  prescribed  d6  /dx  (inverse  mode  for  the  boundary 
layer)  but  also  allows  the  inviscid  flow  to  be  calculated 
with  a  prescribed  geometry  (direct  mode  for  the  inviscid 
flow) . 


5.  Semi-Inverse  Method 

Our  derivation  of  the  semi-inverse  method  differs 
from  that  given  by  Le  Balleur.  Moreover,  as  we  shall  see 
the  formula  we  derive  for  use  in  supersonic  regions  dif¬ 
fers  from  his. 

Let  a®  denote  the  value  of  o  which  is  the  solution 
of  the  viscous-inviscid  coupling  problem.  The  corresponding 
velocity  gradient  will  be  denoted  by  du®/dx.  Then  it  fol¬ 
lows  from  Eq.  (4.7)  that 

a®  =  A  +  B  ^  (5.1) 

Suppose  as  in  our  previous  analysis 

a  =  a®+Aa  =  a^+ee^'^*  (5.2) 


The  corresponding  velocity  gradient  du/dx  calculated  from 
the  boundary  layer  equation  satisfies 

o  =  A  +  B  ^  (5.3) 


Subtract  Eq.  (5.1)  from  Eq.  (5.3)  and  use  Eq.  (5.2).  Then 


ivx  .  e 

ee  ^  ^  _  du 

B  dx  ”  dx 


(5.4) 


For  M<1,  it  follows  from  Eq.  (4.14)  that 

ivx  _  du  du® 

S  dx  dx 


(5.5) 


where  du/dx  is  the  velocity  gradient  calculated  from  the 
inviscid  flow  (at  y  =  0).  Now  subtract  Eq.  (5.5)  from 
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Eq.  (5.4).  Then 


1  u  V  .  j 

ri  _  °°  .  ivx  _  ^  ^ 
''B  B  ^  '  ~  dx  “  dx 


(5.6) 


which  can  be  written  as 


Acx 


ee 


ivx 


BB 

0  -  u„Bv 


rdU 

^dx 


du, 

dx-* 


(5.7) 


If  at  the  n^^  iteration  a  =  a° ,  then  the  (n■^-l)^^ 
iteration  provides  the  required  solution  if  = 

According  to  Eqs.  (5.7)  and  (5.2)  this  will  be  the  case 
if 


n+1 _  n  _  BB  rdu  _  dui 

“  “  u„Bv  -  B  ^dx  dxJ 

00 


(5.8) 


where  M  <  1 . 

Similarly  if  we  start  with  Eq .  (4.18)  and  use  Eq. 
(4.4),  we  find  that  for  M>1 


V  j  ,  .  e 

”  ivx  _  ^  du 

B  dx  ■  dx 


(5.9) 


Subtract  Eq.  (5.9)  from  Eq.  (5.4) 


T  lU  V  .  J 

A  ,  “  ^  IVX  du  du 

^  =  to  -  di 


(5.10) 


In  order  to  eliminate  the  annoying  imaginary  term 
iu^v/B  we  first  differentiate  Eq.  (5.10)  with  respect 
to  x; 


+ 


lu 


■)ee 


IVX 


d^u  _  ^ 

dx^  dx^ 


(5.11) 


Multiply  Eq.  (5.11)  by  u^B/B  and  subtract  from  Eq.  (5.10). 
Then 
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1 

<B*  — 


)Ge 


ivx 


du  _  ^ 
dx  dx 


u  B 
00 


^2~ 

rd_U 

^  dx 


d^u. 

dx  J 


(5.12) 


We  now  use  the  same  reasoning  which  led  to  Eq .  (5.8)  to 
obtain 


n+1  n 
a  -  a  = 


BB 


2„2  2  ^ 

u_B  V  + 


tu^B[ 


dx2 


d  u-i 

2^ 

dx 


_  ft  f^M  _ 


(5.13) 


where  M  >  1 . 
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When  the  error  in  a  is  simple  harmonic  with  frequency 
V,  it  can  be  eliminated,  with  use  of  formulas  (5.8)  and 
(5.13),  for  our  model  problem,  in  one  iteration.  Unfor¬ 
tunately,  however,  the  error  in  a  will  actually  be  a  super¬ 
position  of  harmonics  with  frequencies  v'  limited  only  by 
the  mesh  width  and  the  size  of  the  computational  domain. 
Accordingly  it  behooves  us  to  investigate  the  performance 
of  Eqs.  (5.8)  and  (5.13)  when 


n  _  e  .  „iv ' X 
a  =  a  +  e  e 


(6.1) 
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Since  B<  0  for  Green's  lag-entrainment  method,  it  is 
apparent  that  |p(v')l  <  1  for  all  possible  v'  if  we  croose 
V  to  be  the  maximum  possible  frequency: 


^  =  ^max  = 


With  this  substitution  Eq.  (6.5)  becomes 

$  -  Bu^v  ' 

^  ^  ^  ■  e  -  Bu  V 

“  max 


(6.6) 


(6.7) 


For  attached  flows  it  is  usually  found  that  I Bu  v|  <<  6 
for  all  possible  v  unless  a  very  fine  mesh  is  used.  From 
Eq .  (6.7)  we  see  that  for  attached  flows  we  would  expect 


Vi(v’)  *=  0 


(6.8) 


indicating  that  Eq.  (5.8)  works  very  well  in  this  case. 

For  separated  flows,  on  the  other  hand,  we  can  expect 

1 >>  B  so  that 
'  OO  max' 


Bu^v '  -  B 

1  - 

OO  max 

For  high  values  of  the  frequency  v'  this  becomes 


(6.9) 


y(v')  =  l-v'/v 


max 


while  for  low  frequencies 


m(v')  =  1  + 


Bu  V 
OO  max 


(6.10) 


(6.11) 


(We  recall  that  B < 0  so  that  lp(v')l  <1.)  Hence  in  the 
separated  case  high  frequency  errors  are  reduced  sub¬ 
stantially  with  each  iteration,  while  low  frequency  errors 


are  also  reduced  but  only  slowly.  Experimental  verifica¬ 
tion  of  these  results  will  be  given  later. 
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Now  consider  the  performance  of  the  supersonic  itera¬ 
tion  formula  (5.13)  when  is  given  by  Eq.  (6.1).  It 


follows  from  (5.12)  that 

ufB(v'  )^ 


■)ee 


iv  '  X 


j  U  B  ,2~  .2 

du  _  du  _  00  rd_u  _  d__u, 

dx  ~  dx  8  j  2  j  2-1 

dx  dx 


(6.12) 


Substitute  Eqs. 

n+1 

a 


(6.1)  and  (6.12)  into  Eq. 


+  (1  - 


„2,  2„2.  ,.2 

8  +  u  B  (v  '  ) 

00 

„2  ^  2^2  2 

8  +  u  B  V 


■)e  e 


(5.13). 
iv  '  X 


Then 


(6.13) 


Again  it  follows  that  if  we  choose  v=v  ,  then  errors  of 

max 

all  frequencies  will  be  reduced.  Moreover,  just  as  in  the 
subsonic  case,  high  frequency  errors  are  always  substantially 
reduced,  while  low  frequency  errors  are  reduced  only  slowly 
when  separation  occurs. 

When  actually  implementing  formulas  (5.8)  and  (5.13) 
we  use 


V  =  it/Ax 


u„  =  u 


8  =  ll-M" 


(6.14) 


where  Ax  is  the  local  mesh  width,  u  is  the  local  inviscid 
velocity  and  M  is  the  local  Mach  number.  Formula  (5.8)  can 
be  written  in  a  form  which  agrees  with  that  of  Le  Balleur 
[1980]  if  we  introduce  the  non-dimensional  parameter  B 

* 

B  =  B  (6.15) 

u 


I 


Substitute  Eqs.  (6.14)  and  (6.15)  into  Eq.  (5.8).  Then 
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characteristic  of  separated  flows)  then  Eq.  (6.21)  yields; 

U(v')  =  1  -  (6.22) 


In  this  case  lii(\^')|  <  1  and  we  would  expect  convergence. 
If,  however,  for  some  frequency  V  we  have  | Bu^v ' |  <<  6 
(which  is  characteristic  of  attached  flows),  then  Eq. 
(6.21)  gives 


igBu^v ’ 

y(\)'  )  =  1  +  ”2 - 2 

e  + (u  Bv  ) 

^  '  “  max 


(6.23) 


In  this  case  |u(v')|  >1  so  that  we  would  not  expect 


convergence . 

In  summary  we  would  expect  Le  Balleur's  formula 
Eq.  (6.18)  to  converge  for  separated  or  nearly  separated 
flows,  but  note  that  it  may  fail  to  converge  when  the 
boundary  layer  is  attached. 

Naturally,  both  Eqs.  (5.13)  and  (6.18)  were  tested 
numerically  in  a  transonic  calculation.  In  the  super¬ 
sonic  region  ahead  of  the  shock,  where  the  boundary  layer 
is  firmly  attached,  formula  (5.13)  converged  very  quickly. 
On  the  other  hand,  (6.18)  failed  to  converge.  This 
numerical  test  confirms  the  theoretical  analysis  given 


above . 
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If  we  compare  Eq .  (7.3)  with  Eq.  (5.8)  we  see  that  we  would 
expect  the  Carter  algorithm  to  be  valid  for  subsonic  flows 
if 


<  BB 

u  Bv  -P 
<»  max 


(7.4) 


In  cases  where  Eq.  (7.4)  is  satisfied  strongly  we  would 
expect  the  Carter  algorithm  to  work  even  when  over¬ 
relaxation  is  employed. 

If  we  compare  Eq .  (7.3)  with  Eq.  (5.13)  we  see  that, 


for  supersonic  flow,  the  Carter  algorithm  omits  the  term 
containing  the  expression 


(7.5) 


Accordingly  it  is  of  interest  to  analyse  Eq.  (5.13)  with 
this  term  deleted.  Then 


n+1  n  -BS  rdu  du. 


(7.6) 


If  we  assume  that  a  is  given  by  Eq.  (6.1)  it  follows  from 
Eq.  (5.10)  that 


<h 


^iv '  X  _  ^  ^ 

dx  dx 


(7.7) 


If  we  substitute  Eqs.  (6.1)  and  (7.7)  into  Eq.  (7.6)  we 
find  that  the  original  error  in  a°  is  multiplied  by  the 


factor 


+ iBu^Bv' 

P  ( V  '  )  =  1  ~  n  o  O  O 

8  + 

»  max 


(7.8) 


If  1 u^Bv I «  8  for  all  possible  v  (which  is  characteristic 
of  attached  flows)  then 


p(v ' )  =«  0 


(7.9) 


so  that  the  Carter  algorithm  is  simulating  a  convergent 
algorithm.  On  the  other  hand  if  |  u^^^^Bv  |  >>  S  (which  is 
characteristic  of  separated  flows)  then 


p(v’)  =  1  - 


u  Bv'^ 

®  max 


(7.10) 


In  this  case  |u(v')l  > 1  so  that  we  would  not  expect 
convergence . 

In  summary,  we  would  expect  the  Carter  algorithm  to 
converge  for  attached  supersonic  flows  but  not  for  separated 
supersonic  flows. 
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A  coupling  procedure  employing  Eqs.  (5.8)  and  (5.13) 
is  to  be  preferred  over  the  Carter  algorithm  in  any  case 
because  it  has  a  sounder  theoretical  basis. 
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8.0  Numerical  Tests 


In  order  to  assess  the  performance  of  Eqs.  (5.8)  and 
(5.13),  transonic  calculations  about  an  18%  thick 
circular  arc  airfoil  at  zero  degrees  angle  of  attack  were 
carried  out.  A  Cartesian  computational  grid  of  size 
110 X  30  was  used  for  the  inviscid  analysis.  In  the  x- 
direction  50  equally  spaced  mesh  points  were  placed  on 
the  airfoil,  while  30  mesh  points  were  placed  both  up¬ 
stream  and  downstream  of  the  airfoil,  exponentially 
stretched  to  a  distance  of  five  chord  lengths  away  from 
the  airfoil  in  each  direction.  The  mesh  was  also  expo¬ 
nentially  stretched  in  the  y-direction  to  a  distance  of 
five  chord  lengths  from  the  airfoil.  The  size  of  the  mesh 
in  the  y-direction  near  the  airfoil  was  1%  of  chord.  Free 
stream  conditions  were  assumed  to  prevail  at  the  far  field 
boundaries,  while  the  boundary  condition  at  y = 0  was 
imposed  in  accordance  with  Eq.  (3.42). 

Equations  (5.8)  and  (5.13)  were  implemented  in  their 
integrated  forms.  Given  the  local  inviscid  velocity  u  and 
the  velocity  u  calculated  from  the  boundary  layer,  we 
updated  6  using  the  following  equations.  For  M<  1 


-  (6*)"  = 


u„Bv-8 


[u-u] 


(8.1) 


While  for  M  >  1 
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^  n  +■  1  ^  n 

(6  )  -  (6  )  = 


B3 


2^2  2  2 
u  B  V  +3 


r  ordu  dUi 

'^00  tdx  dx^ 


3  [u  -  u] } 


(8.2) 


We  point  out  for  emphasis  that  v,  and  3  are  given  in 
Eq .  (6.14)  and  B  is  the  coefficient  of  du/dx  occurring 
in  the  boundary  layer  equation.  Green's  lag  entrainment 
method  yields 

B  =  F26*/(Hu)  (8.3) 


The  boundary  layer  was  calculated  in  the  inverse  mode 
(prescribed  d6*/dx)  everywhere  using  the  same  mesh  as 
used  by  the  inviscid  code.  Equation  (2.8)  indicated  that 
for  the  sake  of  numerical  accuracy  in  integrating  the 
boundary  layer  equations,  the  inverse  mode  was  to  be  pre¬ 
ferred  over  the  direct  mode  everywhere  except  in  the  far 
wake  region. 

Once  6*  is  updated  using  Eqs.  (8.1)  and  (8.2)  it  is 

* 

important  to  calculate  d6  /dx  required  by  the  boundary 
layer  equations  in  the  appropriate  manner.  Usually  for 
M  <  1  we  used  a  centered  scheme 


^  =  '^i+1  "  '^i-1 

^  dx  ■'i 


’‘i+l  "  ^i-1 


(8.4) 


However,  at  the  trailing  edge,  6  suffers  a  slope  discon¬ 
tinuity  necessitating  the  use  of  one-sided  differences. 

To  calculate  the  boundary  layer  upstream  of  the  trailing 
edge  (occurring  at  i =  ITE)  we  used 
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(MI) 

''  dx  ''i=lTE 


(8.5) 


While  for  the  wake  calculation  downstream  of  the  trailing 
edge  we  used 


*  ♦ 

,d6  v  _  1+1  1 

^  dx  '^i=ITE  ^i+1  “ 


(8.6) 


The  boundary  layer  equations  were  integrated  up  to  the 

* 

trailing  edge  using  Eq.  (8.5)  to  supply  d6  /dx  at  the 

trailing  edge.  The  boundary  layer  calculation  was  then 

* 

restarted  using  (8.6)  to  give  d6  /dx  at  the  trailing  edge. 
This  procedure  also  allows  us  to  deal  with  the  slope 
discontinuities  due  to  switching  from  the  boundary  layer 
mode  to  the  wake  mode  in  Green's  method. 

For  M  >  1  it  was  found  best  to  calculate  do  /dx  using 
a  difference  scheme  with  upwind  bias.  Consider  the 
formula 


6*  -  6*  6*  -  5* 

dx'i  ■ 


.d6 


“i-1 


1+1 


(8.7) 


Using  X  =  .75  to  introduce  upwind  bias  gave  much  better 
convergence  in  the  supersonic  i:one  upstream  of  the  shock 
than  did  the  centered  scheme  indicated  in  Eq.  (8.4).  On 
the  other  hand,  using  A  =  . 25  which  introduced  a  downwind 
bias  produced  very  poor  results. 

A  series  of  plots  relating  to  test  calculations  appear 
at  the  end  of  the  paper.  The  first  sequence  covers  the 
case  where  =  .7425  and  the  Reynolds  number  based  on  chord 
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length  is  4.10^.  The  standard  Green's  lag  entrainment 

method  was  used  for  this  calculation.  Figure  8.11  shows 

the  pressure  distribution  which  arises  from  a  purely 

inviscid  calculation.  Figure  8.12  shows  the  pressure 

distribution  which  arises  from  viscous-inviscid  coupling. 

Figures  8.11  and  8.12  are  superimposed  in  Fig.  8.13.  As 

can  be  seen,  the  boundary  layer  drives  the  shock  a  small 

distance  upstream  and  slightly  lowers  the  maximum  Mach 

number  which  is  obtained.  Also  it  smooths  out  the  pressure 

* 

distribution  at  the  trailing  edge.  A  plot  of  6  appears 

♦ 

in  Fig.  8.14.  Note  that  6  rises  dramatically  in  response 
to  the  strongly  unfavorable  pressure  gradient  supplied  by 
the  shock  wave.  This  wedging  effect  accounts  for  the 

weakening  of  the  shock  and  its  upstream  displacement. 

* 

Also  note  the  cusp  in  6  which  occurs  at  the  trailing 
edge.  A  graph  of  the  skin  friction  coefficient  appears  in 
Fig.  8.15.  This  drops  dramatically  in  response  to  the 
shock  wave  and  actually  becomes  negative  (indicating 
separated  flow)  near  the  trailing  edge.  Downstream  of  the 
trailing  edge,  Green's  method  is  applied  to  a  wake  calcula¬ 
tion  which  demands  that  the  skin  friction  coefficient  vanish. 
Figure  8.16  compares  the  pressure  calculated  by  the  boundary 
layer  code  with  that  determined  by  the  inviscid  code. 

Ideally  of  course  our  coupling  algorithm  should  ensure  that 
these  are  identical.  However,  as  might  be  expected  some 
trouble  is  experienced  within  the  shock  profile  and  to  a 
minor  extent  near  the  trailing  edge.  Despite  this  dif¬ 
ficulty,  the  pressure  distribution  as  calculated  by  the 
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boundary  layer  is  in  excellent  agreement  with  the  solu¬ 
tions  to  the  Navier-Stokes  equations  presented  by  Deiwert 
[1976] .  Wind  tunnel  interference  corrections  were  not 
applied  to  the  experimental  data.  As  a  consequence,  there 
is  a  small  discrepancy  between  the  calculated  and  observed 
results  (Mehta  and  Lomax  [1981]). 

It  should  be  pointed  out  that  for  this  example  the 
original  Le  Balleur  algorithm  given  by  equation  (6.18) 
fails  to  converge  in  the  supersonic  region  ahead  of  the 
shock.  Le  Balleur  overcomes  this  problem  by  resorting  to 
the  classical  iteration  (suitably  under-relaxed  in 
accordance  with  his  theory  presented  in  section  4)  in 
regions  of  attached  flow.  The  modified  method  on  the 
other  hand,  is  rapidly  convergent  in  this  region  of 
attached  supersonic  flow.  The  decision  as  to  when  one 
should  switch  to  the  classical  iteration  can  now  be  made 
on  the  basis  of  numerical  accuracy  in  computing  the 
boundary  layer  and  need  not  be  influenced  by  considera¬ 
tions  of  coupling  algorithm  convergence. 

As  a  more  severe  test  of  our  coupling  procedures  we 
also  considered  the  case  =  .788.  The  Reynolds  number 

C 

based  on  chord  was  again  taken  to  be  4.10  .  The  standard 
Green's  lag  entrainment  method  was  used.  Figure  8.21  shows 
the  purely  inviscid  pressure  distribution  together  with  the 
pressure  calculated  from  the  viscous-inviscid  coupling. 

In  this  case  the  boundary  layer  has  quite  a  profound  effect 

* 

on  the  outer  inviscid  flow.  The  corresponding  6  and  skin 
friction  coefficient  are  plotted  in  Figs.  8.22  and  8.23. 
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Note  that  we  have  shock  induced  separation,  which  accounts 
for  the  strong  influence  of  the  boundary  layer.  A  con¬ 
vergence  history  of  the  coupling  algorithm  is  presented  in 

V 

Figs.  (8.24)  through  (8.28).  As  predicted  by  the  theory 
presented  in  section  6,  the  algorithm  does  converge  quite 
rapidly  in  regions  of  attached  flow  and  the  rate  of  con¬ 
vergence  is  slower  in  regions  of  separated  flow.  Con¬ 
siderable  difficulty  is  experienced  within  the  shock  pro¬ 
file.  However  the  graceful  recovery  downstream  of  the 
shock  is  a  tribute  to  the  basic  stability  of  the  coupling 
algorithm.  In  order  to  achieve  convergence  for  this  case 
it  was  found  necessary  to  under-relax  the  coupling 
algorithm  in  the  supersonic  region.  A  relaxation  factor 
of  0)  =  .5  was  used,  although  a  larger  value  probably  would 
have  sufficed.  In  the  subsonic  region  a  slight  amount 
of  over-relaxation,  say  <0  =  1.2,  is  generally  helpful, 
although  none  was  employed  in  this  calculation. 

The  pressure  distribution  downstream  of  the  shock 
predicted  by  this  calculation  does  not  agree  well  with 
experiment.  The  experimental  data  indicate  that  there  is 
a  supersonic  plateau  downstream  of  the  shock.  This  example 
was  recalculated  using  the  modified  Green's  lag  entrain¬ 
ment  method  due  to  East  [1977].  The  results  are  presented 
in  Figs.  8.31  through  8.34.  As  can  be  seen,  the  Cp  dis¬ 
tribution  downstream  of  the  shock  is  augmented  and  the 
shock  is  moved  further  upstream.  The  calculated  shock 
location  is  in  excellent  agreement  with  experiment,  but 
the  supersonic  plateau  downstream  of  the  shock  is  still 


not  predicted.  The  Navier-Stokes  solutions  presented  by 
Deiwert  [1976]  experienced  the  same  difficulty.  It  is 
currently  believed  that  the  only  way  to  overcome  this 
problem  is  to  prescribe  the  experimentally  measured 
pressures  at  the  downstream  boundary  (Mehta  and  Lomax 
[1981] ). 

Each  of  the  three  calculations  presented  here  was 
run  for  200  iterations,  although  the  final  answer  was 
obtained,  for  all  practical  purposes,  after  100  itera¬ 
tions.  The  computer  time  required  for  each  case  was 
approximately  18  seconds  on  the  CDC  7600.  Only  about 
20%  of  this  time  was  devoted  to  the  boundary  layer 
calculation.  Thus  the  total  time  required  to  compute  a 
coupled  viscous-inviscid  solution  is  comparable  to  that 
required  to  compute  the  inviscid  flow  alone. 


becomes  large.  Rewrite  Eq .  (8.8)  as 
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^  =  i  _  A 

dx  B  dx  B 


(8.9) 


♦ 

Errors  in  d6  /dx  are  ameliorated  by  the  small  factor  1/B. 

Hence  even  though  the  coupling  algorithm  introduces  errors 
*  ~ 

in  d5  /dx,  the  value  of  du/dx  is  still  close  to  that  which 

* 

would  be  calculated  from  the  exact  value  of  5  .  This 
explains  why  the  pressure  coefficients,  calculated  by  the 
boundary  layer,  always  appear  to  be  so  reasonable,  even 
though  the  inviscid  pressure  coefficients  do  not. 

It  may  be  desirable  to  introduce  more  sophisticated 
procedures  in  the  regions  around  the  shock  and  the  trailing 
edge,  such  as  triple  deck  theory  (see  Inger  [1981]  and 
Melnik  [1980]).  In  this  case  our  coupling  algorithm  pro¬ 
vides  a  stable  framework  in  which  to  embed  these  more 
sophisticated  techniques. 
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9.  Calculations  with  the  Exact  Code 


Encouraged  by  the  success  of  our  numerical  tests  In 
section  8,  we  now  perform  calculations  using  the  exact 
code.  The  effect  of  the  boundary  layer  on  the  Invlscld 
flow  Is  represented  through  the  surface  transpiration 
condition 

^n  =  A  (Pe'^e^*^ 

(see  Eq.  (2.9)).  In  order  to  apply  the  semi-inverse 
method  for  this  case  we  cast  the  boundary  layer  equations 
in  the  form 


(9.2) 


(9.3) 


It  is  then  easily  shown  that  formulae  (5.8)  and  (5.13) 
apply,  provided  that: 


o  is  replaced  by 

B  is  replaced  by 

6  is  replaced  by  Bp_u 

©  © 

As  in  section  8,  the  algorithm  was  implemented  in  its 
integrated  form.  In  Eqs.  (8.1)  and  (8.2)  we  make  the 
replacements  mentioned  above  and  also  note  that 
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6*  is  replaced  by  m 

A  series  of  calculations  were  carried  out  for  an 
NACA  0012  airfoil  at  zero  angle  of  attack.  A  body-fitted 
mesh  of  size  110x30  generated  by  Eisman's  method  [1979] 
was  used.  This  mesh  is  finer  than  the  mesh  displayed  in 
Figs.  3.1  and  3.2.  Eighty  mesh  points  were  placed  on  the 
airfoil.  The  mesh  was  concentrated  near  the  leading  edge 
but  the  largest  mesh  width,  which  occurred  near  the 
trailing  edge,  was  only  1.8%  of  chord.  Thirty  mesh  points 
were  used  in  the  wake,  exponentially  stretched  to  a 
distance  of  five  chord  lengths  away  from  the  airfoil.  The 
mesh  was  also  exponentially  stretched  in  the  pseudo  radial 
direction  to  a  distance  of  five  chord  lengths  away  from 
the  airfoil.  The  size  of  the  mesh  in  the  pseudo  radial 
direction  near  the  airfoil  and  wake  centerline  is  1%  of 
chord.  Free  stream  conditions  were  assumed  to  prevail  at 
the  far  field  boundaries. 

Five  sets  of  calculations  with  corresponding  experi¬ 
mentally  measured  pressure  distributions  are  presented  in 
Figs.  9.11  through  9.55.  The  experiments  were  carried  out 
at  the  National  Aeronautical  Establishment  (NAE)  facility 
in  Ontario,  Canada  (see  Thibert,  Grandjacques  and  Ohman 
[1979]).  The  pressure  coefficients  measured  on  the  lower 
surface  of  the  airfoil  are  plotted  as  O's  while  those  per¬ 
taining  to  the  upper  surface  are  plotted  with  A's.  The 
experiments  were  conducted  at  high  Reynolds  numbers,  so 
in  the  calculations  it  was  assumed  that  the  boundary  layer 
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was  in  a  fully  developed  turbulent  state  at  5%  of  chord. 

The  first  set  of  data  corresponded  to  =  .490.  The 
flow  is  strictly  subsonic  so  that  potential  theory  should 
be  valid.  Indeed  the  purely  inviscid  calculation  is  in 
excellent  agreement  with  experiment.  The  effect  of  the 
boundary  layer  is  totally  negligible  except  very  close 
to  the  trailing  edge. 

The  second  set  of  data  was  run  for  M  =  .693.  The 

oo 

flow  still  remains  subsonic  with  the  maximiim  experimentally 
measured  Mach  nvimber  being  .932.  Thus,  once  again,  po¬ 
tential  theory  should  be  valid  for  the  inviscid  portion  of 
the  flow.  In  this  case,  however,  there  is  a  small  dis¬ 
crepancy  between  calculated  and  measured  pressures  near 
the  leading  edge.  The  effect  of  the  boundary  layer  is 
predicted  to  be  negligible  in  this  region.  Therefore,  the 
coupled  solution  is  in  only  slightly  better  agreement  with 
experiment  than  is  the  purely  inviscid  calculation. 

A  possible  explanation  for  the  discrepancy  between 

theory  and  experiment  is  that  the  boundary  layer  remains 

laminar.  In  the  experiment  the  transition  to  turbulence 

was  left  free  and,  unfortunately,  the  position  of  the  free 

transition  was  not  established.  If  the  boundary  layer 

managed  to  remain  laminar  over,  say,  the  first  20%  of  chord 

then  the  viscous  effects  would  be  more  profound  than 

those  calculated  assuming  the  boundary  layer  to  be  turbulent. 

Evidence  supporting  this  conjecture  is  provided  by  the  third 

set  of  data  run  at  M  =  .776.  In  this  case  the  flow  becomes 

00 


supercritical  and  we  have  a  weak  shock  near  the  leading 
edge.  However,  the  calculated  shock  position  is  10%  of 
chord  aft  of  the  experimentally  measured  position.  More¬ 
over,  the  experiments  show  that  the  pressure  coefficient 

* 

C  just  downstream  of  the  shock  is  nearly  equal  to  C  . 

P  P 

This  indicates  that  viscous  effects  are  seriously  affecting 
the  shock.  The  calculations  based  on  the  turbulent  boundary 
layer  assumption  show  that  the  shock  is  too  weak  to  produce 
a  substantial  wedging  effect  (see  Fig.  9.32).  Therefore, 
a  turbulent  boundary  layer  should  not  cause  a  dramatic 
change  in  the  shock  jump  conditions.  Thus,  we  conclude 
that  the  boundary  layer  may  still  be  laminar  when  it 
interacts  with  the  shock. 

In  the  last  two  data  sets  run  at  M  =  .814  and  = 

.835,  respectively,  the  boundary  layer  is  certainly 
turbulent  by  the  time  it  interacts  with  the  shock. 

However,  the  shock  locations  are  still  not  properly 
predicted  by  the  calculations.  For  the  M^  =  .835  case,  for 
example,  the  purely  inviscid  solution  predicts  the  shock 
location  to  be  at  20%  of  chord  aft  of  that  measured  experi¬ 
mentally  (see  Fig.  9.51).  The  coupled  solution  drives  the 
shock  forward  5%  of  chord  (see  Fig.  9.54)  so  the  difference 
between  theoretical  and  observed  positions  is  still  15%  of 
the  chord  length. 

The  most  likely  cause  of  this  discrepancy  is  explained 
by  Fig.  9.60  taken  from  Lock  [1980].  Here  we  see  that  for 
=  .85,  the  shock  location  predicted  by  the  conservative 
form  of  the  full  potential  equation  does  in  fact  lie  15% 
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of  chord  downstream  of  that  predicted  by  the  Euler  equa¬ 
tions.  Thus  it  appears  that  if  we  are  to  expect  the 
coupled  viscous- inviscid  calculation  to  agree  with  experi¬ 
mental  data  for  all  flows  of  transonic  interest  we  really 
must  solve  the  Euler  equations. 


10.  Conclusion 


A  viscous-inviscid  coupling  procedure,  which  is  a 
modification  of  that  proposed  by  Le  Balleur,  has  been 
presented.  Using  a  model  problem  suggested  by  Le  Balleur 
we  showed  that  the  modified  algorithm  converges  rapidly 
for  attached  flows  and  also  works  well  for  separated 
flows.  This  theoretical  analysis  was  confirmed  by 
numerical  tests  for  a  series  of  challenging  transonic 
calculations  involving  attached  and  separated  flows.  The 
calculations  attest  to  the  remarkable  stability  of  the 
coupling  algorithm,  but  also  point  out  problems  which 
occur  near  the  shock  and  the  trailing  edge.  Nevertheless 
the  method  developed  in  this  thesis  to  compute  coupled 
viscous-inviscid  flows  should  provide  a  suitable  vehicle 
in  which  to  incorporate  more  sophisticated  treatments  of 
these  regions. 
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Fig.  2.1  Van’ af Ion  of  "F^/H  wi+f) 
H  ai  Mach  Numbers  M-0.0, 
1.0,  1.5  for  Siandard  and 
Modified  (3reen‘*5  Laq- 
Enirainmeirt  Method. 
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FIG.  3.1  MESH  GENERATED  ABOUT  NACA  0012  AIRFOIL 
MESH  IS  61  BY  30 
VIEH  OF  ENTIRE  MESH 
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FIG.  8.12  18X  THICK  CIRCULAR  ARC  AIRFOIL 
PRESSURE  AFTER  VISCOUS- I NV I SC ID 
COUPLING  IS  PERFORMED 
HINF- . 7425  RECH-4E06 
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FIG.  8.16  t8X  THICK  CIRCULAR  ARC  AIRFOIL 
VISCOUS  AND  INVISCID  PRESSURE 
AFTER  COUPLING 
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FIG.  8.21  18X  THICK  CIRCULAR  ARC  AIRFOIL  ,  . 

PURE  INVISCID  PRESSURE  DISTRIBUTION 
AND  PRESSURE  AFTER  COUPLING 
HINF>.78d  RECH-«06 
■COUPLED  * INVISCID 
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FIG.  8.25  i8X  THICK  CIRCULAR  ARC  AIRFOIL 
VISCOUS  AND  INVISCID  PRESSURE 
AFTER  SO  ITERATIONS 
NINF-.788  RECH-<4E06 
■VISCOUS  *INVISCID 


FIG.  8.26  I8X  THICK  CIRCULAR  ARC  AIRFOIL 
VISCOUS  AND  INVISCID  PRESSURE 
AFTER  100  ITERATIONS 
NINF-.78d  RECH«^06 
■VISCOUS  *INVISCID 


FIG.  8.27  I8X  THICK  CIRCULAR  ARC  AIRFOIL 
VISCOUS  AND  INVISCID  PRESSURE 
AFTER  158  ITERATIONS 
HINF-.788  RECH««08 
■VISCOUS 


INVISCID 


FIG.  8.28  18X  THICK  CIRCULAR  ARC  AIRFOIL 
VISCOUS  AND  INVISCID  PRESSURE 
AFTER  288  ITERATIONS 
HINF-.788  RECH-<4E86 
■VISCOUS  *INVISCID 


flQ.  0.31  I8S  THICK  CIRCULAR  ARC  AIRFOIL 

PURE  INVISCID  PRESSURE  DISTRIBUTION 
AND  PRESSURE  AFTER  COUPLING  WITH 
EAST'S  MODIFIED  BOUNDARY  LAYER  CODE 
HINF-.7d8  RECH«^E00 


FIG.  8.32  I8X  THICK  CIRCULAR  ARC  AIRFOIL 
OISPUCEMENT  THICKNESS  USING 
HOOIFIEO  BOUNDARY  LAYER  METHOD 
MINF-.788  RECH-^06 


FIG.  8.33  I8X  THICK  CII?CULAR  ARC  AIRFOIL 
SKIN  FRICTION  COEFFICIENT  USING 
NOOIFIEO  BOUNDARY  LAYER  METHOD 
HINF-.78d  RECH-«08 


FIG.  8.34  I8X  THICK  CIRCULAR  ARC  AIRFOIL 
VISCOUS  AND  INVISCID  PRESSURE  AF 
COUPLING  HOOIFIEO  BOUNDARY  UYER 
HINF-.78d  RECH-4E06 
-VISCOUS  -INVISCID 


FIG.  9.M  PURE  INVISCID  PRESSURE 

AND  PRESSURE  AFTER  COUPLING 
H1NF>  ,499  RECH>  I7.5E*G6 
NACA  0012  AIRFOIL 
■COUPLED  MNVISCID 
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FIG.  9.24  PURE  INVISCIO  PRESSURE 

AND  PRESSURE  AFTER  COUPLING 
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FIG.  9.32  DISPLACEMENT  THICKNESS 
MINF-  .776  RECH> 
NACA  0012  AIRFOIL 
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EXPERIMENTAL  AND  PURE  INVISCID  PRESSURES 
NINF-  .814  RECH-  24.7E*06 
NACA  0012  AIRFOIL 
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FIG.  9.52  DISPLACEMENT  THICKNESS 

MINF-  .835  RECH*  24.7E*06 
NACA  0012  AIRFOIL 
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Fid,  $.54  PURE  IKVISCID  PRESSURE 

AND  PRESSURE  AFTER  CODING 
NINF-  .335  RECH>  24.7E*06 
NACA  0012  AIRFOIL,^^,^ 
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Fig.  9.60  NACA  0012  ,  MI NF  =  -65 
A  Com  joarison  hefujeen  Solutions 
of  the  Eul  er  <;»nd  Potential  Flow 
Equations  . 


) 


END 

DATE 

FILMED 

7-81 


DTIC 


